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Abstract 

We give a generalization to an infinite tree geometry of Vidal's infinite time-evolving block 
decimation (iTEBD) algorithm [4 for simulating an infinite line of quantum spins. We numer- 
ically investigate the quantum Ising model in a transverse field on the Bethe lattice using the 
Matrix Product State ansatz. We observe a second order phase transition, with certain key 
differences from the transverse field Ising model on an infinite spin chain. We also investigate 
a transverse field Ising model with a specific longitudinal field. When the transverse field is 
turned off, this model has a highly degenerate ground state as opposed to the pure Ising model 
whose ground state is only doubly degenerate. 

1 Introduction 

The matrix product state (MPS) description |9j has brought a new way of approaching many-body 
quantum systems. Several methods of investigating spin systems have been developed recently 
combining state of the art many-body techniques such as White's Density Matrix Renormalization 
Group [6] [7] (DMRG) with quantum information motivated insights. Vidal's Time Evolving Block 
Decimation (TEBD) algorithm [TJ [2] uses MPS and emphasizes entanglement (as measured by the 
Schmidt number), directing the computational resources into that bottleneck of the simulation. It 
provides the ability to simulate time evolution and it was shown that MPS-inspired methods handle 
periodic boundary conditions well in one dimension [10 , areas where the previous use of DMRG 
was limited. TEBD has been recast into the language of DMRG in [8 j and adapted to finite systems 
with tree geometry in (3j . DMRG is especially successful in describing the properties of quantum 
spin chains, the application of basic DMRG-like methods is limited for quantum systems with 
higher dimensional geometry. New methods like PEPS |12j generalize MPS to higher dimensions, 
opening ways to numerically investigate systems that were previously inaccessible. 

We are interested in investigating infinite translationally invariant systems. Several numerical 
methods to investigate these were developed recently. The iTEBD algorithm jl] (see also Secj4]) is a 
generalization of TEBD to infinite one-dimensional systems. A combination of PEPS with iTEBD 
called iPEPS [13] provides a possibility of investigating infinite translationally invariant systems 
in higher dimensions. Our contribution is a method to investigate the ground state properties 
of infinite translationally invariant quantum systems on the Bethe lattice using imaginary time 
evolution with Matrix Product States. The Bethe lattice is an infinite tree with each node having 
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Figure 1: The Bethe lattice (infinite Cayley tree). 



three neighbors, as depicted in FigjT] It is translationally invariant in that it looks the same at 
every vertex. This geometry is interesting, because of the following connection to large random 
graphs with fixed valence. Moving out from any vertex in such a random graph, you need to go a 
distance of order logn, where n is the number of vertices in the graph, before you detect that you 
are not on the Bethe lattice, that is, before you see a loop. 

We choose to investigate the quantum transverse field Ising model on the Bethe lattice. Note 
that we work directly on the infinite system, never taking a limit. First we test the iTEBD method 
on a system with a known exact solution, the infinite line. Then we turn to the Bethe lattice with 
the new method we provide. In both cases, the Hamiltonian is given by 

^ = ^E( 1 -^) + 2E( 1 -<')' W 

(m) * 

where the sum over i is over all sites, and the sum over is over all bonds (nearest neighbors). 
We show that imaginary time evolution within the MPS ansatz provides a very good approximation 
for the exact ground state on an infinite line, giving us nearly correct critical exponents for the 
magnetization and correlation length as we approach the phase transition. We obtain new results 
for the quantum Ising model in transverse field on the infinite tree. Similarly to the infinite line, 
we observe a second order phase transition and obtain the critical exponent for the magnetization, 
Pt ~ 0.41 (different than the mean-field result). However, the correlation length does not diverge 
at the phase transition for this system and we conjecture that it has the value l/ln2. 

We also investigate a model where besides an antiferromagnetic interaction of spins we add a 
specific longitudinal field \a\ for each spin: 

#notoo = J^(l + < + a! + ( rX)+^(l-4). (2) 

(hi) i 

We choose the longitudinal field in such a way that the interaction term in the computational basis 
takes a simple form, 1 00) (00|j-, giving an energy penalty to the 1 00) state of neighboring spins. 
(We follow the usual convention that spin up in the z-direction is called 0.) We call it the NOT 
00 model accordingly. This model is interesting from a computational viewpoint. The degeneracy 
of the ground state of -ff N oT00 at /i = is high for both infinite line and infinite tree geometry of 
interactions. We are interested in how our numerical method deals with this case, as opposed to 
the double degeneracy of the ground state of ([TJ at h = 0. We do not see a phase transition in this 
system as we vary J and h. 

The paper is organized as follows. Section [2] is a review of the MPS ansatz and contains its 
generalization to the tree geometry. In Section [3j we review the numerical procedure for unitary 
updates and give a recipe for applying imaginary time evolution within the MPS ansatz. In Section 
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[4j we adapt Vidal's iTEBD method for simulating translationally invariant one-dimensional systems 
to systems with tree geometry. Section[5]contains our numerical results for the quantum Ising model 
in a transverse field for translationally invariant systems. In Section 5.1 we test our method for 
the infinite line, and in Section [5 .2| we present new results for the infinite tree. We turn to the not 
00 model in Section [6] and show that our numerics work well for this system even when there is 
a high ground state degeneracy. In Section [7] we investigate the stability of our tree results and 
conjecture that they may be good approximations to a local description far from the boundary of 
a large finite tree system. 



2 Matrix Product States 

If one's goal is to numerically investigate a system governed by a local Hamiltonian, it is convenient 
to find a local description and update rules for the system. A Matrix Product State description is 
particularly suited to spin systems for which the connections do not form any loops. Given a state 
of this system, we will first show how to obtain its MPS description, and then how to utilize this 
description in a numerical method for obtaining the time evolution and approximating the ground 
state (using imaginary time evolution). We begin with matrix product states on a line (a spin 
chain) , and then generalize the description to a tree geometry. In [3j we give a numerical method 
of updating the MPS description for both real and imaginary time simulations. 

2.1 MPS for a spin chain 

Given a state \ip) of a chain of n spins 

W) = Yl c -,s l , Si+ i,..\s 1 ) 1 ...\si) i \s i+1 ) l+1 ...\s n ) n , (3) 

...SiS i + 1 ... 

we wish to rewrite the coefficients c si Sn as a matrix product (see [5] for a review of MPS) 

_ V \(i-l)p(j): s i \ (*) r (*+l),s;+i ( A\ 

t-...,Sj,s i+ i,... — ' " ' a a,b A b 1 b,c A c ••• \^) 

...abc... 

using n tensors and n — 1 vectors A". The range of the indices a,b, . . . will be addressed later. 
After decomposing the chain into two subsystems, one can rewrite the state of the whole system 
in terms of orthonormal bases of the subsystems. \® is the vector of Schmidt coefficients for the 
decomposition of the state of the chain onto the subsystems 1 . . . i and i + 1 . . . n. 

In order to obtain the A's and the T's for a given state one has to perform the following 
steps. First, perform the Schmidt decomposition of the chain between sites i—1 and i as 

l^) = El^)w-l A ^ ) l^)v..,n' ( 5 ) 
a=l 

where the states on the left and on the right of the division form orthonormal bases required to 
describe the respective subsystems of the state \ip). The number Xi-i (the Schmidt number) is the 
minimum number of terms required in this decomposition. The Schmidt decomposition for a split 
between sites i and i + 1 gives 

l^>=El^>W^ ^>*4.i,....»- ( 6 ) 
6=1 
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Figure 2: Two successive Schmidt decompositions on a line allow us to find the T tensor for the 
marked site and the two A vectors for the bonds coming out of it. 



These two decompositions (see FIGj2J) describe the same state, allowing us to combine them to 
express the basis of the subsystem i, . . . ,n using the spin at site i and the basis of the subsystem 
i + 1 , . . . , n as 

\^i,.,n = EE r S SA ?l^l^ + l,..,n> (7) 
8=0,1 6=1 

where we inserted the Al for convenience. This gives us the tensor T®. It carries an index s 
corresponding to the state |s) of the i-th spin, and indices a and b, corresponding to the two 
consecutive divisions of the system (see FIGj2|. Because \4> a ) (and are orthonormal states, 
the vectors A and tensors T obey the following normalization conditions. From ^ we have 

f>P = l, (8) 

6=1 

while ([7]) implies 

<<M0a>v..,n = E ErgrA?r«' s A« = 8 a , a , , (9) 

s=0,l 6=1 

and 

(en^w = E E ^t^Sr^t 1 ^ = v ■ (io) 

s=0,l 0=1 



2.2 MPS on Trees 

Matrix Product States are natural not just on chains, but also on trees, because these can also 
be split into two subsystems by cutting a single bond, allowing for the Schmidt-decomposition 
interpretation as described in the previous section. The Matrix Product State description of a 
state of a spin system on a tree, i.e. such that the bonds do not form loops, is a generalization of 
the above procedure. Tree-tensor-network descriptions such as ours have been previously described 

mm- 

(k) 

Specifically, for the Bethe lattice with 3 neighbors per spin, we introduce a vector Xa k for each 
bond k and a four-index (one for spin, three for bonds) tensor Tal]a t ,a m for each site i. We can then 
rewrite the state \ip) analogously to ([3]),(|4]) as 

w = ( II E A £ } )( II E T S:t,a n )\---)\si)\...), (ii) 

\ fcebonds <Jfc=l / \ iesites Sj / 
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Figure 3: The three Schmidt decompositions on a tree required to obtain the T tensor for the 
marked site and the three A vectors for the bonds emanating from it. 



where ai,a m , a n are indices corresponding to the three bonds I, m and n coming out of site i. Each 
index ai appears in two T tensors and one A vector. To obtain this description, one needs to perform 
a Schmidt decomposition across each bond. This produces the vectors A™. To obtain the tensor 

for site i, one needs to combine the three decompositions corresponding to the bonds of site i 
as depicted in Figj3} Analogously to expressing the orthonormal basis for the first subsystem 
marked in Figj3]in terms of the state of the spin |sj) and the orthonormal bases for the latter two 
subsystems in PigjsJ one obtains the tensor Tai,a m ,a n for site i. 

The normalization conditions for a MPS description of a state on a tree are analogous to Q- ( 10 ) . 
We have 

£ A if = 1 > (12) 



Xk Xl 



VVVrW' s * \«2\(02 r «,s = x /ion 

/ j / j a k ,ai,a m , /y a k /v a ; 1 a k ,ai,a m u a m ,a m , ; V iL V 

s=0,l a k =l a; = l 



and two other variations of (13) with k,l and m interchanged 



3 Simulating Quantum Systems with MPS 

We choose to first describe the numerical procedures for a chain of spins. Then, at the end of the 
respective subsections, we note how to generalize these to tree geometry. 



3.1 Unitary Update Rules 

The strength of the MPS description of the state lies in the efficient application of local unitary 
update rules such as U = e~ lA ^ 1 (where A is an operator acting only on a few qubits). First, 
we describe the numerical procedure in some detail, and then, in the next Section, discuss how to 
modify the procedure to also implement imaginary time evolution. 

Given a state as a Matrix Product State, we want to know what happens after an application 
of a local unitary. In particular, for a 1-local U acting on the i-th spin, it suffices to update the 
local tensor 

r (i),s U TT s r (i),s' / -i j \ 

L a,b >U s' L a,b ■ V i4 J 

The update rule for an application of a 2-local unitary V acting on neighboring spins i and i + 1, 
requires several steps. First, using a larger tensor 

e& = E (rS^M;: 1 ^) a^), (15) 

b 
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we rewrite the state \ip) as 



After the application of V, the tensor in the description of 1^) changes as 

©^E^X/" (17) 



One now needs to decompose the updated tensor to obtain the updated tensors r^, r^ +1 ^ and 
the vector A®. We use the indices a,s and c, t of to introduce combined indices (as) and (ct) 
and form a matrix T^ as \^ with dimensions 2%j_i x 2\i+i as 

Using the singular value decomposition (SVD), this matrix can be decomposed into T = QAW, 
where Q and W are unitary and A is a diagonal matrix. In terms of matrix elements, this reads 

T(as),(ct) = ^ Q (as), b D b,bWb, (ct) ■ (19) 
b 

The diagonal matrix D = diag(A^) gives us the updated Schmidt vector A^*'. The updated tensors 
and r( l+1 ) can be obtained from the matrices Q, W and the definition of (15) using the old 



vectors A*-* ^ and A^ +1 ^ which do not change with the application of the local unitary V. After 



these update procedures, the conditions (|8j)-( 10 ) are maintained. 

The usefulness / succintness of this description depends crucially on the amount of entanglement 
across the bipartite divisions of the system as measured by the Schmidt numbers Xi- To exactly 
describe a general quantum state of a chain of n spins, the Schmidt number for the split through 
the middle of the chain is necessarily Xn/2 = 2 n / 2 . Suppose we start our numerical simulation in 
a state that is exactly described by a MPS with only low Xi s - The update step described above 
involves an interaction of two sites, and thus could generate more entanglement across the i, i + 1 

(i) 

division. After the update, the index bin X b would need to run from 1 to 2\i to keep the description 
exact (unless \i already is at its maximum required value Xi = 2 min l* ,n ~ *1). This makes the number 
of parameters in the MPS description grow exponentially with the number of update steps. 

So far, this description and update rules have been exact. Let us now make the description an 
approximate one (use a block-decimation step) instead. First, introduce the parameter %, which 
is the maximum number of Schmidt terms we keep after each update step. If the amount of 
entanglement in the system is low, the Schmidt coefficients A^ decrease rapidly with b (We always 
take the elements of A sorted in decreasing order). A MPS ansatz with restricted Xi = X wm 
hopefully be a good approximation to the exact state However, we also need to keep the 

restricted x throughout the simulation. After a two-local unitary update step, the vector A® can 

(i) 

have 2x entries. However, if the b > x entries in A^ after the update are small, we are justified 
to truncate A® to have only x entries and multiply it by a number so that it satisfies Q. We 
also truncate the V tensors so that they keep dimensions 2 x x x X- The normalization condition 
(10) for rw will be still satisfied exactly, while the error in the normalization condition ([9l will be 



small. This normalization error can be corrected as discussed in the next section. This procedure 
keeps us within the MPS ansatz with restricted x- 
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The procedure described above allows us to efficiently approximately implement local unitary 
evolution. To simulate time evolution 



\m) 



-iHt 



1-0(0)) 



(20) 



with a local Hamiltonian like ((!]), we first divide the time t into small slices At and split the 



Hamiltonian into two groups of commuting terms i?l and Hm' ■ Each time evolution step e 
can then be implemented as a product of local unitaries using the second order Trotter-Suzuki 
formula 



r(«) 



-WAt 



u 2 



-iH 



(x) At 




-iff 



(x) At 



(21) 



The application of the product of the local unitaries within each group can be done almost in 
parallel (in two steps, as described in Section [4]), as they commute with each other. 



These update rules allow us to efficiently approximately simulate the real time evolution (20) 
with a local Hamiltonian H for a state \ip) within the MPS ansatz with parameter x- The number 
of parameters in this MPS description with restricted x 1S then n(2% 2 ) for the tensors and 
(n — l)x for the vectors A®. The simulation cost of each local update step scales like 0(x 3 ), 
coming from the SVD decomposition of the matrix O. For a system of n spins, we thus need to 
store 0(2nx 2 + nx) numbers and each update will take 0(nx 3 ) steps. 

The update procedure generalizes to tree geometry by taking the tensors F with dimensions 



2 x X x X x X as rn Section 2.2 For a local update (on two neighboring spins i and i + 1 with bonds 
labeled by l,m,n and n,o,p) we rewrite the state analogously to (16) as 



Z-^t Z-^ (a k ai),(a a p ) 
a k ,ai,a ,a p s,t 



J a kl 



K) \ s )i \t)i+l \4>a D ) \(pa p ) 



using the tensor 
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(a k ai),(a a p ) 



\(k) X (l)ST> ( r (A),s A (m) r (i+1),< \\(o)\(p) 
/x a k /K ai / j I 1 a k ,ai,a m "a m a m ,a ,a p I /x a a /x a p > 



(22) 



(23) 



with combined indices (ai-ai) and {a a p ). One then needs to update the tensor O as described 
above (17)-(19). The decomposition procedure to get the updated vector A^ m ^ and the new tensors 
rW and r( l+1 ^ now requires 0(x 6 ) computational steps. The cost of a simulation on n spins thus 
scales like 0(nx 6 )- 



3.2 Imaginary Time Evolution 

Using the MPS ansatz, we can also use imaginary time evolution with e~ Ht instead of (20) to 



look for the ground state of systems governed by local Hamiltonians. One needs to replace each 
unitary term e _Jj4A * in the Trotter expansion (21) of the time evolution with e~ A/ ^ f followed by a 



normalization procedure. However, the usual normalization procedure for imaginary time evolution 
(multiplying the state by a number to keep = 1) is now not enough to satisfy the MPS 



normalization conditions (|8|)-( 10 ) for the tensors T and vectors A we use to describe the state 



The unitarity of the real time evolution automatically implied that the normalization conditions 
(|8j) ,( 10 ) were satisfied after an exact unitary update. While there already was an error in ^ 
introduced by the truncation of the x + 1 • • • %X entries in the non-unitarity of imaginary time 
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evolution update steps introduces further normalization errors. It is thus important to properly 
normalize the state after every application of terms like e~ A ^ t to keep it within the MPS ansatz. 

In [3] , Vidal dealt with this problem by taking progressively shorter and shorter steps At during 
the imaginary time evolution. This procedure results in a properly normalized state only at the end 
of the evolution, after the time step decreases to zero (and not necessarily during the evolution). 



We propose a different scheme in which we follow each local update e 



-A At 



by a normalization 



procedure (based on Vidal's observation) to bring the state back to the MPS ansatz at all times. 
The simulation we run (evolution for time t) thus consists of many short time step updates e~ HAt , 
each of which is implemented using a Trotter expansion as a product of local updates e~ A ^ 1 . Each 
of these local updates is followed by our normalization procedure. 

We now describe the iterative normalization procedure in detail for the case of an infinite 
chain, where it can be applied efficiently, as the description of the state requires only two 
different tensors T (see Section 4.1). One needs to apply the following steps over and over, until 
the normalization conditions are met with chosen accuracy. 

First, for each nearest neighbor pair i, i + 1 with even i, combine the MPS description of these 
two spins (15)-(16), forming the matrix T (18). Do a SVD decomposition of T (19) to obtain a 
new vector AW. The decomposition does not increase the number of nonzero elements of AW , as 
the rank of the 2\ X %X matrix T ( 19 ) was only x (coming from ( 15 )) . We thus take only the first x 

1. using this new AW, we obtain tensors 



(02 



values of AW and rescale the vector to obey Yla=i *a 
r w , r( l+1 ) from (19), and truncate them to have dimensions x x X x 2. Second, we repeat the 
previous steps for all nearest neighbor pairs of spins i, i + 1 with i odd. 

We observe that repeating the above steps over and over results in exponential decrease in 
the error in the normalization of the T tensors. We note though, that the rate of decrease in 
normalization errors becomes much slower near the phase transition for the transverse field Ising 
model on an infinite line (see Section 5.1 ). 

In practice, we apply this normalization procedure by using the same subroutine for the local 
updates e _AA *, except that we skip the step (17), which is equivalent to applying the local update 
with At = 0. The normalization procedure is thus equivalent to evolving the state repeatedly 
with zero time step (composing two tensors T and decomposing them again) and imposing the 
normalization condition on the vectors A. Note though, following from the definition of the SVD, 
that each decomposition assures us that one of the conditions (|9]),(10) is retained exactly for the 
updated tensors T. The errors in the other normalization condition for the T tensors are decreased 
in each iteration step. 

The numerical update rules for a system with tree geometry are a simple analogue of the update 



(A) fi- 
nales for MPS on spin chains. Every interaction couples two sites, with tensors T a and T 



c,d,e ' 



with the three lower indices corresponding to the bonds emanating from the sites. One only needs 



to reshape the tensors into r|^' s c and 



(B),t 

c,(de) 



and proceed as described in ( 15 ) and below. 



4 MPS and Trans lationally Invariant Systems 
4.1 An Infinite Line 

For systems with translational symmetry such as an infinite line all the sites are equivalent. We 
assume that the ground state is translationally invariant, and furthermore pick the tensors r w 
and vectors AW to be site independent. For fixed x t ne number of complex parameters in the 
translationally invariant MPS ansatz on the infinite line scales as 2% 2 . 
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Figure 4: The parametrization and update rules for the infinite line. 



When using imaginary time evolution to look for the ground state of this system, within this 
ansatz, it is technically hard to keep the translational symmetry and the normalization conditions 
after each update. Numerical instabilities plagued our efforts to impose the symmetry in the 
procedures described above. In |4j, Vidal devised a method to deal with this problem. Let us break 
the translational symmetry of the ansatz by labeling the sites A and B as in FIGJ4] This doubles 
the number of parameters in the ansatz. The state update now proceeds in two steps. Let the site 
pairs AB interact and update the tensors T^ A \ and the vector \( AB \ Then let the neighbor 
pairs BA interact, after which we update the tensors T^ B \ T^ A ' and the vector \( BA \ What we 
observe is that after many state updates the elements of the resulting and the tensors 
differ at a level which is way below our numerical accuracy (governed by the normalization errors) 
and we are indeed obtaining a translationally invariant description of the system. 

One of the systems easily investigated with this method (iTEBD) is the Ising model in a 



transverse field (37) on an infinite line. Vidal's numerical results for the real time evolution and 
imaginary time evolution |4. of this system show remarkable agreement with the exact solution. 
We take a step further and also numerically obtain the critical exponents for this system. Further 
details can be found in Section |5j where we compare these results for the infinite line to the results 
we obtain for the Ising model in transverse field on the Bethe lattice. 

4.2 An Infinite Tree 

For the infinite Bethe lattice, our approach is a modification of the above procedure introduced by 
Vidal. In order to avoid the numerical instabilities associated with imposing site-independent T 
and A after the update steps, we break the translational symmetry by labeling the "layers" of the 
tree A and B (denoted by half-circles and triangles), as in FIGJ5J 

The Bethe lattice is also symmetric under the permutation of directions. Tensors T with full 
directional symmetry obey r aifejC = T b ^ a = T c ^ b = r CifejQ = r 6jaiC = r a>C)6 . However, for the purpose 
of simple organization of interactions, we will also partially break this symmetry by consistently 
labeling an 'inward' bond for each node, as denoted by the flat sides of the semi-circles and the 
longer edges of the triangles in Fig|5j This makes the first of the three indices of T a fi )C special. 
However, we keep the residual symmetry T a ^ yC = T a ^ b . This we can enforce by interacting a 
spin with both of the spins from the next layer at the same time. The update procedure for the 
interaction between the spins now splits into two steps, interacting the layers in the AB order first, 
and then in the BA order as in FigJHJ Similarly to what we discovered for the line, the differences 
in the elements of the final and are well below the numerical accuracy of our procedure. 

The scaling of this procedure is more demanding than the 0(x 3 ) simulation for a line. The 
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Figure 5: The two-layer, directed labeling of the tree. 





Figure 6: The two-step interactions for the infinite tree. 



number of entries in the matrix used in each update step is 2% 2 x 4% 4 , therefore the SVD 
decomposition requires 0(x 8 ) steps. The scaling of our numerical method is thus 0(x S ) f° r each 
update step. 

4.3 Expectation Values 

A nice property of the MPS state description is that it allows efficient computation of expectation 
values of local operators. First, for a translationally invariant system on a line (with only one 
tensor T and one vector A), we have for an operator acting only on the i-th spin 



E 



x x 



(24) 



=0,1 



o=l 6=1 
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where , = (s^|0®|sj). Similarly, for the expectation values of O^O^ (assuming j > i), 



T YYo il) ,o ij) , 



(25) 



s i ,s^...,s J ,^. a,e b,b>. 



x(^r:: 6 A b r^A c --.A d rJ e A e ). 



Defining a x 2 x X 2 matrix B (where one should think of (bb') as one combined index ranging from 
1 to x 2 ) as 



B{bb'),{cc') - y^r^ e rg* e ,A c A c /, 

s 

and vectors v and w with elements again denoted by a combined index (bb') = 1 ... x 2 as 



(dd') 



J J 



ej j 



we can rewrite ( 25 ) as 



v T BB---Bw, 

j-i-l 



(26) 

(27) 
(28) 

(29) 



There is a relationship between the eigenvalues of the matrix B and the correlation function 
(O^O^) — (OW)(OW). One of the eigenvalues of B is \i\ = 1, with the corresponding right 
eigenvector 



fifed) - X]y^ r b,c r b*c( A c) 2 ' 



and left eigenvector 



(30) 



(31) 



which can be verified using the normalization conditions ([9]) and (10). We numerically observe that 
\i\ = 1 is also the largest eigenvalue. (Note that > 1 would result in correlations unphysically 
growing with distance.) Denote the second largest eigenvalue of B as [i<i- Using the eigenvectors of 
B, we can express B^ 1 ^ 1 in (29), as 

B j-i-l = {lL)p{XR)T + ^-lp { 2L )f3 (2R)T + (32) 

When computing the correlation function, the term that gets subtracted exactly cancels the leading 
term involving /xi = 1. Therefore, if I//2I is less than 1, (32) implies 



(33) 



The correlation function necessarily falls of exponentially in this case, and the correlation length £ 
is related to /i2 as £ = — l/ln/X2- 
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The computation of expectation values for a MPS state on a system with a tree geometry can 
be again done efficiently For single-site operators the formula is an analogue of (24) with 
three A vectors for each T tensors which now have three lower indices. For two-site operators, the 
terms in (29) now become 



B(cc'),(dd') 



v (cc>) 



W(dA>\ 



E^ ) S 'E r 2e, / rK, / (A e ) 2 (A /) ^ 



(34) 
(35) 

(36) 



The correlation length is again related to the second eigenvalue of the B matrix as in ( 33 ) . 



5 Quantum Transverse Field Ising Model 

Our goal is to investigate the phase transition for the Ising model in transverse magnetic field ([I]) 
on the infinite line and on the Bethe lattice. We choose to parametrize the Hamiltonian as 

^=|E( 1 -^i) + ^T^E( 1 -4), (37) 

(hi) * 
where < s < 1 and 6 is the number of bonds for each site (6 = 2 for the line, 6 = 3 for the tree). 



We will investigate the ground state properties of (37) as we vary s. The point s = corresponds 
to a spin system in transverse magnetic field, while s = 1 corresponds to a purely ferromagnetic 
interaction between the spins. 



5.1 The Infinite Line. 

We present the results for the case of an infinite line and compare them to exact results obtained via 
fermionization (see e.g. |14j . Ch.4). Vidal has shown |4 that imaginary time evolution within the 
MPS ansatz is capable of providing a very accurate approximation for the ground state energy and 
correlation function. We show that even using x smaller than used in jl], we obtain the essential 
information about the nature of the phase transition in the infinite one-dimensional system. We 
also obtain the critical exponents for the magnetization and the correlation length. 

In FIGjTJ we show the how the ground state energy obtained using imaginary time evolution 
with MPS converges to the exact energy as \ increases. The exact solution for a line has a second 
order phase transition at the critical value of s, sl = §, and the ground state energy and its first 
derivative are continuous, while the second derivative diverges at s = sl- We plot the first and 
second derivative of E with respect to s obtained numerically and compare them to the exact values 
in FIG. [8j observing the expected behavior already for low \- 

The derivative of the exact magnetization M = (a z ) is discontinuous at sl, with the magneti- 
zation starting to rise steeply from zero as 

Moz(x-x L f, (38) 
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Figure 7: Transverse Ising model on an infinite line. Fractional difference of the ground state energy 
obtained using MPS and the exact ground state, near the phase transition at sl = §• The energy 
scale is logarithmic. 
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Figure 8: Transverse Ising model on an infinite line. The first and second derivative with respect 
to s of the ground state energy obtained via MPS compared with the exact result. 
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Figure 9: Transverse Ising model on an infinite line. Magnetization obtained using MPS vs s, with 
B z = 10" 8 . 



with the critical exponent 0i = 
the transverse field strength in (37) and is 



Here x is the ratio of the ferromagnetic interaction strength to 



2(1 



(39) 



with the value x = xl = 1 at the phase transition (sl = |). We plot the magnetization obtained 
with our method in FIGJ9J To obtain the magnetization depicted in the plot, we used a small 
symmetry breaking longitudinal field with mag nitude B z = 1CT 8 . In FIGmH we plot M vs. x — xl 
on a log-log scale. We also plot a line with slope 0.125. Observe that as x increases, the data 
is better represented by a line down to smaller values of x — xl- For the largest x we display, a 
straight line fit of the data between 4 x 10~ 3 < x — xi < 10" 1 gives us a slope of 0.120. Note that 
the mean field value of the critical exponent for magnetization is 0.5. 



The correlation function (<rfVf ) 



can be computed efficiently using (|25|). Away 



from criticality, it falls off exponentially as e -1 * - ^'^ The falloff of the correlation function is 
necessarily exponential as long as fi2, the second eigenvalue of B, is less than 1. The exact solution 
for the correlation length £ near the critical point has the form 



XL 



(40) 



as x approaches xl- Already at low x the iTEBD method captures the divergence of the correlation 
length. In FIG. 11, we plot £ vs. \x — xl\ on a log-log plot, together with a line with slope — 1. 



Again, as x increases, the data is better represented by a line closer to the phase transition. For 



the highest x we display, a straight line fit of the data between 2x1 < xl — x < 4 x 10 
gives us a slope of —0.92. Note that the mean-field value of the critical exponent for the correlation 
length is 0.5 (corresponding to slope —0.5 in the graph). 
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Figure 10: Transverse Ising model on an infinite line. Log-log plot of magnetization vs. x — xl- 
We also plot a line with slope (3l = 0.125. 
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Figure 11: Transverse Ising model on an infinite line. Log-log plot of the correlation length vs. 
\x — xl\. We also plot a line with slope — = —1. For each \ i n the plot, we choose a numerical 
value of the critical point sl as the point at which the correlation length is maximal. 
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Figure 12: Transverse Ising model on an infinite tree. The first and second derivative with respect 
to s of the ground state energy obtained via MPS. 



5.2 The Infinite Tree (Bethe Lattice). 



The computational cost of the tree simulation is more expensive with growing x than the line 
simulation, so we give our results only up to \ = 8. We run the imaginary time evolution with 
10000 iterations (each iteration followed by several normalization steps) for each point s, taking 
a lower x result as the starting point for the procedure. We also add a small symmetry-breaking 
longitudinal field with magnitude B z = 10~ 8 . 

We see that the energy and its first derivative with respect to s are continuous. However, we 



now observe a finite discontinuity in the second derivative of the ground state energy (see FIG. 12 ), 
as opposed to the divergence on the infinite line. This happens near s = st ~ 0.5733. 

Similarly to the one-dimensional case, the magnetization quickly grows for s > st, while it has 
(nearly) zero value for s < st (see FIG. 13). In FIG. 14 we plot the magnetization vs. x — xt on 
a log- log scale, where x is 



3(1-*)' 



with the value x = xt ~ 0.451 at the phase transition (where s 
whether the magnetization behaves like 



M oc (x — xt) 



13 



(41) 



st ~ 0.5733). We want to test 



(42) 



for x close to xt, which would appear as a line on the log-log plot. As x grows, the data is 
better represented by a straight line closer to the phase transition. If we fit the x = 8 data for 
4 x 10~ 4 < x — xl^x < 10~ 3 , we get f3 = 0.41. We add a line with this slope to our plot. Note 
that the mean-field value for the exponent (3 is 0.5, just as it is for the infinite line. 
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Figure 13: Transverse Ising model on an infinite tree. Magnetization vs. s, with B z 
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Figure 14: Transverse Ising model on an infinite tree. Log-log plot of magnetization 
with B z = 10 -8 . We also plot a line with slope 0.41. 
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Figure 15: Transverse Ising model on an infinite tree. A linear plot of the correlation length vs. s. 



We observe that the correlation length now rises up only to a finite value (see FIG. 15). As 
we increase %, the second eigenvalue of the B matrix, /j,2, approaches a maximum value close to \. 
We conjecture that the limiting value of [12 is indeed I, which corresponds to a finite correlation 
length with value (ln2) _1 . Note that for the infinite line, the second eigenvalue of B approaches 1, 
and so the correlation length is seen to diverge at the phase transition. 



6 The Not 00 Model 

We now look at a model with a different interaction term. Starting with an antiferromagnetic inter- 
action, we add a specific longitudinal field at each site. As in the previous section, we parametrize 
our Hamiltonian ^ with a single parameter s: 

iWoo = *£|(i+*i+^+*ta)+^^£(i-<4), ( 43 ) 

(trfV v ' i 

with b = 2 on the line and b = 3 on the tree. We choose the longitudinal field in such a way that the 
nearest-neighbor interaction term Hj becomes a projector, expressed in the computational basis 
as 

Hij = |00) (00|y , (44) 
thus penalizing only the 1 00) configuration of neighboring spins. Accordingly, we call this model 



NOT 00. The ground state of the transverse Ising model (37) at s = 1 has degeneracy 2. For (43) 
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Figure 16: The not 00 model on an infinite line. The ground state energy and its first two 
derivatives with respect to s. 



on the infinite line or the Bethe lattice, the degeneracy of the ground state at s = 1 is infinite, as 
any state that does not have two neighboring spins in state |0) has zero energy. 

6.1 Infinite Line 



We use our numerics to investigate the properties of (43) on the infinite line as a function of s 



Our numerical results show continuous first and second derivatives of the energy with respect to 



s (see FIG. 16). The magnetization M = (a z ) decreases continuously and monotonically from at 



s = to a final value of —0.606 at s = 1 (see FIG. 17). The second eigenvalue of the B matrix (33) 



rises continuously from at s = 0, approaching 0.603 at s = 1 (see FIG 17). Because ^ < 1> the 
correlation length £ is finite for all values of s in this case. These results imply that there is no 
phase transition for this model as we vary s. 

As a test of our results, we compute the magnetization at s = 1 exactly for this model on a 
finite chain (and ring) of up to n = 17 spins. We maximize the expectation value of Hb = °f 



X 



within the subspace of all allowed states at s = 1 (with no two zeros on neighboring spins), thus 



minimizing the expectation value of the second term in (43) for s approaching 1. We compute the 
magnetization M = {a\) for the middle i = |_^J spin for the ground state of the NOT 00 model 
exactly for a finite chain and ring of up to 17 spins at s = 1. As we increase n, the value of M 
converges to —0.603 (much faster for the ring, as the values of M for n = 14, 17 differ by less than 
10 -4 ). Recall that we obtained M = —0.606 from our MPS numerics with x = 16 for the NOT 
00 model on an infinite line. We also compare the values of the Schmidt coefficients across the 
central division of the finite chain {n = 16) to the elements of the A vector obtained using our MPS 
numerics with \ = 32. We observe very good agreement for the 11 largest values of A^, with the 
difference that our MPS values keep decreasing (exponentially), while the finite-chain values flatten 
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Figure 17: The NOT 00 model on an infinite line. Magnetization as a function of s and correlation 
length as a function of s. 



out at around Xk>u ~ 10 9 (see FIG .18). The behavior of the components of A from MPS doesn't 
change with increasing \- 



In the ground state of (43) at s = 1, the overlap with the 1 00) state of any two neighboring 
spins is exactly 0. If the T tensors are the same at every site, the component of the state \tp) that 
has overlap with the state 1 00) on nearest neighbors can be expressed as 

^(A a r° 6 A b r° c A c )l^)|OO)|0 c ). (45) 

a,b,c 

Furthermore, when the T tensors are symmetric, the elements of the A vectors must be allowed 
to take negative values to make this expression equal to zero. Note that until now, we used only 
positive A vectors, knowing that they come from Schmidt decompositions, which give us the freedom 
to choose the components of A to be positive and decreasing. 

The negative signs in the A vector can be absorbed into every other T tensor, resulting in a 
state with two different T (for the even and odd-numbered sites) and only positive A's. In fact, this 
is what we observe in our numerics, which assume positive A, but allow two different T tensors (see 



4.1). If we allow the elements of A to take negative values, our numerically obtained T tensors are 



identical. 



6.2 Infinite Tree 



Here, we numerically investigate the not 00 model (43) on the Bethe lattice. As on the line, the 



numerics show continuous first and second derivatives of the energy with respect to s (see FIG. 19) 



and a continuous decrease in the magnetization from at s = to —0.671 at s = 1 (see FIG, 20) 



The correlation length behaves similarly as on the line, increasing with s, but it reaches a maximum 



at s = 0.96 for x = 8. The maximum value of £ is apparently lower than l/ln2, (see FIG. 20) 
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Figure 18: Ground state of the NOT 00 model on a line at s = 1. Comparison of the exact Schmidt 
coefficients for a division across the middle of a finite chain and of the MPS values (x = 32) for an 
infinite line. 
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Figure 19: The not 00 model on an infinite tree. The ground state energy and its first two 
derivatives with respect to s. 
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Figure 20: The not 00 model on an infinite tree. Magnetization as a function of s and correlation 
length as a function of s. 



meaning that on the tree, the correlation function (a z a 
faster than 2 - ' 1- - 7 ' for all s. 



{oz )(oz ) falls off with distance 



7 Stability and Correlation Lengths on the Bethe Lattice 

We have found that on the Bethe lattice, for both our models the second eigenvalue y,i of the matrix 
B (33), which determines the correlation length, apparently is never greater than h. In this section 



we argue that this is a model- independent, and calculation method independent, consequence of 
assuming that a translation-invariant ground state is the stable limit of a sequence of ground states 
of finite Cayley trees as the size of the tree grows. For a related problem, the stability of recursions 
for Valence Bond States on Cayley trees has been investigated by Fannes et.al. in 

The Hamiltonians ([I]) and ^ each consist of sums of terms Hj^ , flm , as i n (21 ), where each 



term depends on a single a x and each Hm' on a neighboring pair of a z . We calculate the 
quantum partition function 



r(«) 



Z{0) =tre 



-PH 



(46) 



as the limit of 



Z(N, At) = tr 



-AtH, 



N 



(47) 



as At — > 0, N — > oo with NAt = (3. To find the properties of the ground state, we take (5 — > oo 
so that we need Z(N, At) as At 0, N -> oo with NAt oo and N(At) 3 (to make the 



error in using the Trotter-Suzuki formula go to zero). We interpret (47) as giving the classical 
partition function of a system of Ising spins (s = ±1) on a lattice consisting of N horizontal 
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Figure 21: A system of classical spins on a lattice whose layers are Cayley trees. A and B denote 
the Boltzmann factors. 



layers, each of which is a Cayley tree of radius M (i.e with a central node and concentric rings of 
3, 3 x 2, 3 x 2 2 , . . . , 3 x 2 M - X nodes). We can write d47} as 



Z(N,At) = J2U A U B ^ ( 48 ) 

W 

where the sum is over all configurations of N x (3 x 2 M - 2) spins s = ±1 and the products are of a 
Boltzmann factor A for each horizontal link in the Cayley trees, and a Boltzmann factor B for each 



vertical link between corresponding nodes in neighboring layers (see FIG. 21) (layer N is linked to 
layer 1 to give the trace). The factor A for the link between nodes i,j in the same horizontal layer 
is given by 

A( Si , Sj ) = e - AtH $) {si ' Sj) . (49) 
The factor B for the link between nodes i, i' in the same vertical column is given by 

B( 8i ,ati = (a z = S ;|e~ A * H « \a z = Si ). (50) 

(z) (x) 
When all the terms H). i are of the same form, as are all the terms H).J , the form of the factors A 

and B does not depend on which particular links they belong to. 

Each term in the sum, divided by Z, can be thought of as the probability of a configuration 

{s}. In what follows we will keep N and At fixed and consider the limit M — * oo, i.e. finite Cayley 

tree — * Bethe lattice. We will then suppose that our results, which are independent of the form of 

A and B (provided A, B > 0) will also hold after the iV — > oo limit is taken, i.e. for the quantum 

ground state. 

We will think of our lattice as a single tree, with each node being a vertical column of iV spins. 
(A recent use of this technique to investigate the quantum spin glass on the Bethe lattice is in [16J.) 
Denote by s the vector of N values of s along a column. Let K{s) be the product of the N factors 
B(s, s') along a column and let L(s, s') be the product of ./V factors A(s, s') on the horizontal links 
between nearest neighbor columns. Let Zm be the partition function for a tree of radius M. We 
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Figure 22: The M = 2 Cayley tree. 



can calculate Zm by a recursion on M as follows: 

Z M = ^2K(s)[F M (s)]\ (51) 

s 

F M (s) = ^2l(s,s')K(s')[F M -i(s')] 2 , (52) 

s' 

F (s) = 1. (53) 



It is easy to see that this recursion gives the correct Zm (the case M = 2 is shown in FIG 22 ). We 
also see that 

Pm(s) = -^K(s)[F M (s)] 3 (54) 

is the probability of the configuration s along the central column. For the case N = 1, i.e. classical 
statistical mechanics on a tree, this is the well-known method to find an exact solution |15j . 

In order to have a well-defined translationally invariant limit as M — > oo we would like the 



recursion (52) for Fm to have an attractive fixed point F which Fm approaches as M —> oo. 
'Attractive' means that if we start the recursion with a different Fq(s), sufficiently close to Fq(s) = 
1, the limiting value of Fm(s) will be the same fixed point. This in turn implies that on a Cayley 
tree with large M, small changes in the Hamiltonian on the outer edge will have small effects on 
the properties of the central region. 



First however we need to fix the overall normalization of Fm{s), since if Fm(s) satisfies (52), 
so does a 2 Fm{s) which rules out an attractive fixed point. 
Let 

F M (s) = Z M F M (s), (55) 

so that 

Y,K(s)[Fm(s)) 3 = 1, (56) 
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and 

P M {s) = K{s)[F M {s))\ 

The recursion relation becomes 

hi(s) = A M £) L (^0*(0[i^-i(0] 2 . 



with Am determined by the normalization condition (56). We can now suppose that 



with 



and 



Fm(s) -> F(s) as M -> oo, 

= A^L(s , s')i<'(s / )[F(s / )] 2 , 

^K( S -)[F( S -)] 3 = 1. 



To determine whether F is an attractive fixed point, let 

F M (s) = F( s -) + / A/ (s), 
Am = A(1 + 6m), 



with /m — > and — > as M — > oo. To first order in /m, £Mi (58) and (|56j) become 

f M (s) = e M F(s)+2Y,T(s,s')fM~i(s'), 

s' 

J2k(s)IHs)} 2 Im(s) = 0, 



where 



From (60), 



T(s,s') = XL(s,s')K(s')F(s'). 



^T(s,s')F(s') = F(s), 



and since L(s, s') = L(s', s), 



s)}\ 



(57) 
(58) 

(59) 
(60) 

(61) 



(62) 
(63) 



(64) 
(65) 

(66) 
(67) 

(68) 



i.e. the linear operator T has an eigenvalue one, with right eigenvector F and left eigenvector KF 2 

(69) 



(which from ( |61[ ) have scalar product one). (64) now gives 

J2W)[Hs)] 2 fM(s) = e M + 2j2K(s)[Hs)?fM-i(s), 
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so from (65), €m = 0. Let 



T^s, s') = T(s, s') - F(s)K(s') [F(s')} 2 



(70) 



so that 



Ys^is^s'Ms') = 0, 



(71) 



and 



(72) 



(64) now becomes 



f M (s) = 2^t ± (s,^')/m-i(s')- 



(73) 



( 73 ) shows that F(s) is an attractive fixed point if and only if 

ll^ll < ^ (74) 

(for a tree with valence p+ 1 at each vertex, | is replaced by -). 

We can in fact prove that there does exist an F(s) satisfying (60) and (61), for which the 
corresponding T 1 - has a maximum eigenvalue less than Define a function of F(s) by 



m = ^[F(s)] 2 K(s)L(s,s')K(s')[F(s'tf 



(75) 



We look for a maximum of $ with F(s) restricted to the region 

"£k(s)[F(s)} 3 = 1, 

s 

F(s) > 0. 



(76) 
(77) 



A maximum must exist, but it might be on the boundary of the region, i.e. it might have F(s) = 
for some values of s . Elementary calculations (omitted here) establish that stationary values of $ 
on the boundary cannot be maxima. At stationary points in the interior of the region, i.e. with 
F(s) > for all s, (60) and (61) must be satisfied. If such a stationary point is a maximum, all 
the eigenvalues of 1 — 2T 1 - are > 0, i.e. all the eigenvalues of T 1 - are < g- This is weaker than the 
attractive fixed point condition, which also requires that no eigenvalue is less than — ^, but does 
correspond to our observed property of fi2- 

We now examine the joint probability distribution of sq and Sd, where denotes the central 
column and d a column distance d from the center. From FIGJ23I we see that 



PM(so,s d ) 



1 

Zm 



[Fm{s )} 2 K{s )L(s , si)F m _ 1 (s 1 ) . . . 
• • • L(s d -x, s d )K(s d )[F M „ d (s d )} 2 . 



(78) 



SI,— > s d-l 
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Figure 23: Computing the probability distribution of sq and s d . 



As M — » oo (with <i fixed) this becomes (using (66)) 

P(s ,s d ) = C[F(s )] 2 K(s )T d (s ,s d )F(s d ), 
where the normalization C is determined by 

^P(s ,s d ) = l. 



Using (67) we find 



P(s ) = J2P(s ,s d ) = CK(s )[F(s )] 3 



(79) 



(80) 



(81) 



so from (61), C = 1 (and P(sq) agrees with the limit of (|57[)). Expressing (79) in terms of T 1 - 
using ((70), (Pfll and ((721), 



P(s , r d ) - P(r )P(5d) = K(so)[F(so)] 2 (T ± ) d (so, s d )F(s d ). 



(82) 



Thus the correlation between s*o and s d falls off as /U rf , where \i is the eigenvalue of T 1 - with maximum 
modulus, and so from (74), faster than l/2 rf . 



If this conclusion is correct (and clearly the argument is less than rigorous), it establishes more 
than our experimental observation that \X2 < \- The quantum limit of P(so,s d ) encodes not only 
the static correlation {tpo\ sos d \ipo) in the ground state, but also the imaginary time dependent 
correlation (^g| e Ht soe~ Ht s d \ipo) which in turn determines the linear response as measured by s d 
to a time-dependent perturbation proportional to So- If this indeed falls off faster than l/2 d , then 
we can have some hope that the Bethe lattice can be used as a starting point for investigation of 
fixed valence random lattices. 
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